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ABSTRACT 

We observe two secondary eclipses of the strongly irradiated transiting planet WASP-33b, in the K s 
band at 2.15 /im, and one secondary eclipse each at 3.6 /im and 4.5 /zm using Warm Spitzer. This planet 
orbits an A5V <5-Scuti star that is known to exhibit low amplitude non-radial p-mode oscillations at 
about 0.1% semi-amplitude. We detect stellar oscillations in all of our infrared eclipse data, and also in 
one night of observations at J-band (1.25 /im) out of eclipse. The oscillation amplitude, in all infrared 
bands except K s , is about the same as in the optical. However, the stellar oscillations in K s band 
(2.15 /im) have about twice the amplitude (0.2%) as seen in the optical, possibly because the Brackett-7 
line falls in this bandpass. As regards the exopla netary eclipse , we us e our best-fit values for the eclipse 



depth, as well as the 0.9 /im eclipse observed by ISmith et al. (12011), to exp l ore p ossible states of the 
exoplanetary atmosphere, based on the method of iMadhusudhan fc Seagen (|2009f ). On this basis we 



find two possible states for the atmospheric structure of WASP-33b. One possibility is a non-inverted 
temperature structure in spite of the strong irradiance, but this model requires an enhanced carbon 
abundance (C/O > 1). The alternative model has solar composition, but an inverted temperature 
structure. Spectroscopy of the planet at secondary eclipse, using a spectral resolution that can resolve 
the water vapor band structure, should be able to break the degeneracy between these very different 
possible states of the exoplanetary atmosphere. However, both of those model atmospheres absorb 
nearly all of t he stellar irradiance w ith minimal longitudinal rc-distribution of energy, strengthening the 
hypothesis of ICowan &; Agol (|2011h that the most strongly irradiated planets circulate energy poorly. 
Our measurement of the central phase of the eclipse yields e cosuj = 0.0003 ± 0.00013, which we regard 
as being consistent with a circular orbit. 

Subject headings: stars: planetary systems - transits - techniques: photometric 



1. INTRODUCTION 

Extrasolar planets that orbit close to their stars are 
subject to an intense flux of stellar irradiation. The ro- 
tation of a very close-in planet is expected to become 
tidally locked to its orbital period o n an astrophysically 
short time scale (|Guillot et al.lll996[) . Consequently, the 
most close-in exoplanets will receive stellar irradiation ex- 
clusively on their star-facing hemispheres. The resulting 
heating is believed to b e distributed by strong zonal winds 
(Sho wman et alj |2008). but the dynamics of the zonal re- 
distribution, and therefore the overall energy budget of 
the planet, are affected by the vertical temperature struc- 
ture of the planetary atmosphere. Many close-in planets 
exhibit inverted temperature structures (|Knutson et al.l 
120081 ISeage? k, Demiiigll2010h . probably driven by radia- 
tive absorption in a h igh altitude layer of the atmosphere 
( Bur rows et al.l [20071 ) . The nature of the absorber has 



been actively discussed (jFortnev et al J 120081 ; ISpiegel et al.l 
2009), but remains unknown. 

One promising avenue of investigation is to look for 
correlations between the planetary temperature inver- 
sion and the stellar flux at ultraviolet (UV) wavelengths 
(jKnutson et al.ll2010l ) . Stellar UV radiation has the poten- 
tial to dissociate absorbing molecular species, and to create 
(or d estroy) absorbers via photochemistry (Za hnle et al.l 
|2009[ ) . The spectral distribution of UV flux may be criti- 
cal to the inversion phenomenon. Therefore it is desirable 
to investigate planets orbiting strong sources of far-UV ra- 
diation (i.e., magnetically active stars), as well as planets 
receiving irradiation by thermal UV radiation (i.e., hot 
stars). 

An important planet in the latter category is WASP- 
33b, that orbits an A- typ e J-Scuti star with an or- 
bital period of 1.22 days (jCollier-Cameron et al.l 120101 ; 
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iHerrero et al.l 1201 it) . The large radius and high tem- 
perature of an A-type star produce stronger irradiance 
than would be the case for a solar-type star. Although 
only an upper limit is available f or the mass of WASP- 
33b (|Coilier-Cameron et al.l [2010l ). the planet is impor- 
tant because it is among the most strongly irradiated 
planets. In contrast to other strongly irradiated plan- 
ets such as WASP - 12. (ICowan et al.ll2Q12l: [Crossfield et al 



2012t IZhao et al.l 120121: ICampo et all 1201 It I Croll et al 



201l dMadhusudhan et al.ll2011a ). little is currently known 



about the response of WASP -33b to the strong stellar ir- 
radiation. ISmith et al.l f|201lh measured the thermal emis- 
sion of WASP-33b from secondary eclipse observations at 
0.9 /jm, but there are currently no reported detections of 
the planet at wavelengths longward of 1 /im. 

Stellar intensity oscillations of WASP-33 are seen 
with 0.1% semi-a mplitude at optical wavelengths 
(jHerrero et al.ll20TTl ). The stellar oscillations may exhibit 
a greater or lesser amplitude at infrared (IR) wavelengths. 
The dependence of the oscillation amplitude on wave- 
length potentially carries information on the physics of 
the oscillations in the stellar atmosphere. 

In this paper, we report measurement of the thermal 
emission from WASP-33b, based on ground-based obser- 
vations of two secondary eclipses in the K s band (2.15 /um), 
space-borne observations of eclipses at 3.6- and 4.5 /im by 
Warm Spitzer, and measurement of the intensity oscilla- 
tions of the star in all these IR bands, as well as in J-band 
(1.25 /xm). We describe the observations and extraction 
of photometry from the data in Sec. 2 & 3. In Sec. 4, 
we analyze the data to determine the parameters of the 
planet's eclipse and the oscillatory properties of the star. 
We explore and discuss the implications of our results in 
Sec. 5, and Sec. 6 summarizes our results. 

2. OBSERVATIONS 

2.1. Ground-based Observations 

We observed secondary eclipses of WASP-33b on UT 10 
and 16 October, 2011, using the FLAMINGOS infrared 
HgCdTe imager at the Cassegrain focus of the 2.1-meter 
telescope at Kitt Peak National Observatory. The sky on 
both nights was cloudless, with excellent photometric con- 
ditions prevailing, especially on 10 October. The obser- 
vations used a K s filter, and we defocused the telescope 
so that the diameter of stellar images - measured at the 
half-intensity level - was 30-pixels (18 arc-sec). This sub- 
stantial defocus improves the photometric precison and 
photon-collection efficiency for this bright star (WASP-33 
has V=8.3). We compensated for drift in the telescope de- 
focus using manual updates at approximately 30 minute 
intervals, based on a known formula that relates focus po- 
sition to temperature and zenith distance. In that way, we 
attempted to maintain the greatest possible image stabil- 
ity. 

The 20 x 20 arc-min field of view of FLAMINGOS 
(2048 x 2048 pixels) provided 5 comparison stars imaged 
simultaneously with WASP-33. We obtained a nearly con- 
tinuous sequence of 20-second exposures on each night, 
amounting to 725 exposures on 10 October and 574 expo- 
sures on 16 October. Including the overhead of reading the 
detector and writing FITS files, the observational cadence 
was 45 seconds per exposure. We verified that the times 



written to the FITS headers were free of clock errors (to 
about 1 second precision) by comparing the header values 
to manual timing made using a web-displayed UT clock. 
In addition to WASP-33, we observed multiple sky expo- 
sures with position offsets, that were used to construct 
a flat-field frame by median-combining the offset sky ex- 
posures. We also observed dark frames using the same 
exposure times as for the stellar and sky-flat observations. 

All of our WASP-33 observations used a quad-detector 
off-axis guide camera, sensitive to optical radiation, and 
producing real-time pointing corrections for the telescope. 
The guider was very effective at damping image motion on 
time scales of minutes, but differential refraction between 
the optical and the infrared leads to a slow drift in stellar 
positions, amounting to about 1 arc-sec over 60-degrees of 
zenith distance. This slow positional drift has only a small 
effect on our photometry. 

During the observations, we noted a significant instabil- 
ity in the infrared signal from the FLAMINGOS detector, 
that occurred in response to the changing position of the 
telescope. In order to diagnose and correct for these insta- 
bilites, we obtained a sequence of K-band exposures during 
the day, with the dome closed and dark, and the primary 
mirror cover closed. These 2-second exposures measured 
the thermal emission and scattered light from the tele- 
scope mirror cover, and we moved the telescope to differ- 
ent positions because the signal instabilities seemed to be 
a function of telescope position, as described in Sec. 3.2. 

During our WASP-33 oservations, we also noticed rel- 
atively prominent intensity oscillations of the star in K s - 
band. Hence, we also observed 651 consecutives exposures 
of WASP-33 on the night of UT 14 October, when no 
eclipse (or transit) of the planet occurred. This time se- 
ries was observed using J-band, and helps to establish the 
degree to which the amplitude of stellar oscillations in in- 
tensity may depend on wavelength. 

2.2. Warm Spitzer Observations 

Warm Spitzer observed one eclipse of WASP-33 at 
3.6 /im, for 9.6 hours on 26 March 2011, as well as obser- 
vations of equal duration spanning one eclipse at 4.5 /im 
on 30 March 2011. The observations were made under the 
Cycle-6 Target of Opportunity program (J. Harrington, 
P.I.). Both observations used subarray mode to collect 
1264 data cubes, each data cube comprising 64 exposures 
of a 32 x 32-pixel section of the detector. The exposure 
time was 0.4 seconds per exposure, and the observations 
were not interrupted by any re-pointing. 

3. PHOTOMETRY 

3.1. Ground-Based Photometry 

In both K s and J-band, we subtract a dark frame from 
each image, and divide by the sky-flat for that wavelength. 
Each sky flat is normalized to have an average value near 
unity; this facilitates conversion of the observed images 
from data numbers to electrons using the known gain of 
the detector electronics (4.9 e/DN). In the limit of large de- 
focus, the stars are essentially images of the telescope's pri- 
mary mirror; they are not well-approximated using Gaus- 
sians or other functions commonly used for centroiding. 
Therefore we determine the center position of each star by 
exploiting the sharp edges that are characteristic of the 
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pupil images. Summing the image of each star in the X- 
coordinate produces intensity as a function of Y, I(Y). 
The sharp edges of the pupil image will produce peaks in 
the derivative dl/dY, one peak at each edge of the image. 
We find those peaks in the spatial derivative, and adopt 
an average of the Y-positions of derivative peaks on each 
side of the I(Y) profile as being the center of the star in Y, 
and vice- versa for X. With the center of each star deter- 
mined for each exposure, we calculate the flux in a circular 
aperture of a given radius centered on that star. 

In addition to WASP-33, we measure 4- to 6 comparison 
stars on each image. The identification of the comparison 
stars are given in Table 1, with their J and K magnitudes. 
The subset of comparison stars that were actually used 
varied from night to night, due to using the J-band wave- 
length on 14 October, and K-band sky transparency and 
thermal background that varied between 10 and 16 Octo- 
ber. The subset of comparison stars used for each night 
are noted in the caption of Table 1. 

For both the target and comparison stars, we deter- 
mined the value of the sky background by constructing a 
histogram of pixel values in a 100 x 100-pixel box surround- 
ing each star. The sky background dominates the peak of 
those histograms, and we fit a Gaussian to each histogram, 
and thereby determine the sky background value for each 
star as the centroid of each Gaussian. Multiplying the 
background value per-pixel for each star, times the area of 
the aperture containing that star, yields the background 
contribution for that star. Subtracting the background 
from each aperture measurement yields the stellar intensi- 
ties of WASP-33 and the comparison stars for that image. 

We repeat our photometry by varying the radius of the 
synthetic aperture from 18 to 30 pixels, in 1-pixel incre- 
ments. The 13 different aperture radii produce 13 different 
realizations of photometry for each night. Our rationale 
for varying the aperture radius is to optimize the signal- 
to-noise and the match between the target star and the 
comparison stars. Although our images are strongly defo- 
cused, there is still scattered light beyond the edges of the 
defocused pupil image for each star. Variable seeing and 
errors in aperture centering cause the total flux intercepted 
by a chosen aperture to vary. Measurements with larger 
apertures are less sensitive to these variations, but include 
more background noise. The optimum aperture radius is 
approximately 20 pixels, but we determine it separately 
for each night, as described in Sec. 4.1. 

3.2. Instrumental Instabilities 

As noted in the Introduction, we experienced instabil- 
ities in the photometric signals. These instabilities were 
discovered using quick-look evaluation of the data during 
each night of observations. The nature of the instabilities 
is that sharp increases or decreases in stellar intensities 
occur typically 4 to 6 times each night. These changes 
in intensity are large compared to our photometric pre- 
cision: sudden changes as large as 4% were seen. Stars 
close together in angular extent experienced similar, but 
not identical effects. For our observations, the comparison 
stars were typically several hundred pixels distant from 
WASP-33, so the instabilities were not precisely common- 
mode. The sky background also exhibited this instability, 
but to a much smaller relative degree - barely detectable. 



During our observing run, we noticed that the signal in- 
stabilities never occurred for stars at negative declination, 
and they tended not to occur at large hour angles. We con- 
cluded that these puzzling instabilities were triggered by 
motion of the telescope, and that motivated our diagnostic 
observations described in Sec. 2. The root cause of the sig- 
nal instabilities is now known: following our observing run, 
Dick Joyce of the Kitt Peak scientific staff disassembled 
FLAMINGOS, and discovered that a 4-40 screw (proba- 
bly from the filter box) had come loose and fallen onto the 
first camera lens, producing field-dependent vignetting as 
it shifted position. We were initially tempted to discard 
all of these 'screwy' data, but the excellent photometric 
quality of the nights (especially on 10 October) motivated 
us to develop a methodology to remove the effect of the 
rolling screw. 

Fortunately, the nature of this effect - a sudden shift 
in position of the screw followed by periods of stability - 
makes it amenable to robust correction using only the data 
themselves. We used our test observations with the dome 
closed to validate our correction methodology. We produce 
aperture photometry from these frames using apertures of 
20-pixel radius, and no background subtraction (the back- 
ground is the signal for the test observations). We center 
the apertures at the same locations as WASP-33 and the 
comparison stars, and extract time series signals. Figure 1 
(top panel) shows the time series at the WASP-33 posi- 
tion. The large and sudden changes in signal level are 
obvious, and they correspond to times of telescope mo- 
tion. The telescope motion for these test observations was 
done in a series of 0.5-hour increments, versus continuous 
tracking for the actual stellar observations. Nevertheless, 
the results are very similar, and we conclude that we have 
successfully reproduced the signal instabilities. 

We correct these instabilities by operating on the time 
derivative of the signal. We calculate the distribution of 
the ensemble of the numerical time derivative values for 
each time series. This distribution is primarily Gaussian, 
with extreme outliers that correspond to the times of sig- 
nal instability. We fit to the Gaussian core of the dis- 
tributions, to determine the unbiased standard deviation 
of the signal derivative at the position of each star. We 
adopt a threshold value (typically 3<r) and we correct the 
outlying points in each series of derivative values, where 
the absolute value of the derivative exceeds this threshold. 
This threshold value (in a) is an adjustable parameter in 
our correction procedure. After identifying the derivative 
points beyond the threshold, one option is to set these out- 
lying derivative points to zero. However, in practice we 
replace derivative outliers with a 5-minute smoothed ver- 
sion of the derivative. We choose the 5-minute smoothing 
time because it represents the typical time for significant 
changes in stellar intensity, e.g., as caused by the telluric 
atmosphere. After correcting the outlying points in the 
derivatives, we integrate the corrected derivatives to yield 
corrected time series. 

Figure 1 (upper panel) shows the corrected time series 
for our dome test observations at the WASP-33 position. 
The derivative of the time series is shown in the middle 
panel, with the rejected outliers marked. The lower panel 
of Figure 1 shows the result of dividing the corrected time 
series at the WASP-33 position, by the total of the cor- 
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rected series at the comparison star positions. In the ab- 
sence of our correction procedure this quotient would show 
considerable noise at the several-percent level. The cor- 
rection procedure successfully produces a quotient whose 
short-term variations are of order 200 parts-per-million 
(ppm), with long-term variations approaching 700 ppm. 
We do not regard the long-term variations as meaningful 
to our correction procedure, because the background in 
the telescope dome is not guaranteed to be stable at this 
level. We conclude that our procedure to correct for the 
signal instabilites is effective, and we apply it to our stellar 
photometry. 

3.3. Warm Spitzer Photometry 

Photometry of Warm Spitzer data is now a fa- 
miliar exercise dHebrard et al.l l2010t iBeerer et al.l 1201 It 



Deming ct al. 2011a; Dc morv et al 



120111: iDesert et al.l 



20111 iTodorov et al.ll2012D . and we used well-tested pro- 
cedures applied to the BCD files produced by version 
S18.18.0 of the Spitzer pipeline. We performed aperture 
photometry on each frame, in each data cube, for both 
eclipses. We found the centroid of the stellar image by 
fitting a 2-D Gaussian, and computed the flux in a circu- 
lar aperture centered on the star, for aperture radii from 
2.0 to 5.0 pixels in 0.5-pixel increments. We accounted for 
the background contribution in each frame; the per-pixel 
background was measured as the centroid of a Gaussian 
fit to a histogram of pixel values for that frame. Multiply- 
ing the background per-pixel times the area of the aper- 
ture yielded the background contribution that was sub- 
tracted from the flux in the aperture. Aperture radii near 
3.0-pixels produced the highest signal-to-noise photome- 
try, as judged by the point-to-point scatter. Accordingly, 
we adopted a 3.0-pixel radius for our photometry at both 
3.6- and 4.5 /im. 

4. ANALYSIS OF THE PHOTOMETRY 
4.1. Ground-based Photometric Analysis 
4.1.1. Optimization of Photometric Parameters 

As described in Sec. 3.1, we generated 13 versions of 
our ground-based photometry using a range of aperture 
radii, and we also have a free parameter (the threshold) to 
correct the photometry for instrumental instabilities. We 
determine which combination of aperture radius and cor- 
rection threshold produces the most robust results. We 
increment the aperture radius in 1-pixel steps, and the 
correction threshold in steps of 0.1a. For each combina- 
tion of aperture radius and correction threshold, we cal- 
culate the linear Pearson correlation coefficients between 
the photometric time series for WASP-33, and the cor- 
responding photometry for each comparison star. In the 
limit of very high correction threshold (i.e., no instabil- 
ity correction) , there is very little correlation between the 
comparison stars and WASP-33. That is consistent with 
the behavior of the instrumental instability, as evaluated in 
our closed-dome test observations, wherein we found that 
different portions of the detector exhibit different instabil- 
ities. Similarly, we expect photometric aperture radii that 
are too small or too large would degrade the correlation 
between WASP-33 and the comparison stars. 

In order to select the best combination of aperture ra- 
dius and correction threshold, we average the correlation 



coefficients that we compute for the comparison stars ver- 
sus WASP-33. For our J-band observations on 14 Octo- 
ber, the highest average correlation coefficient (0.85) is 
achieved when the combination of (radius, threshold) is 
(22 pixels, 3.1er). For fC s -band observations on 10 and 16 
October, the optimum combinations are (22, 2.8cr) and 
(20, 3.1cr), producing average correlation coefficients of 
0.88 in both cases. These results are quite reasonable, 
because aperture radii near 20-pixels are modestly greater 
than the radius to the half-intensity point of the defocused 
image (about 15 pixels), and 20-pixels is what our intuition 
told us to choose for quick-look analyses at the telescope. 
Similarly, threshold values near 3er are reasonable because 
they allow sudden spikes in the signal derivatives to be 
identified without perturbing normal fluctuations due to 
photometric noise. We conclude that the resultant time se- 
ries photometry is the best that can be achieved from these 
data. As a by-product of this process, we evaluate the 
point-to-point fluctuations in the WASP-33 photometry, 
from the fit to the Gaussian distribution of signal deriva- 
tive values. These fluctuations are sub-milli-magnitude for 
all three nights; we find J-band noise of 642 ppm on 14 Oc- 
tober, and if s -band noise of 725 and 905 ppm on 10 and 
16 October, respectively. The higher noise on 16 Octo- 
ber is produced by a higher thermal background, due to 
warmer weather and slightly degraded (but still excellent) 
atmospheric transparency over Kitt Peak. 

4.1.2. Normalization using Comparison Stars 

Following selection of the best photometric aperture ra- 
dius, and threshold for instrumental correction, we fur- 
ther correct the WASP-33 time series using the comparison 
stars. Normally this would be accomplished by dividing 
WASP-33 by the sum of the comparison stars. However, 
we find improved results using a slightly different proce- 
dure: instead of a straight sum of the comparison stars, we 
use a weighted sum. We denote the intensity of WASP-33 
at time index i as Wi, and we write: 



JV 



(1) 



where is the intensity of the j-th comparison star at 
time index i, and the ctj coefficients - one for each of the 
TV comparison stars - are determined by linear regression 
(matrix inversion). The linear regression seeks to produce 
the best overall match between the weighted sum of the 
comparison stars and WASP-33, i.e., to make Eq. (1) be 
exact. Because of noise and intrinsic variations in WASP- 
33, Eq.(l) can never be solved exactly, only for the set of 
otj that produces the best approximation. Having solved 
for those a.j, we divide the Wi by the right hand side of 
(1). This division by the weighted sum of the comparison 
stars removes instrumental and telluric effects, but leaves 
noise and the intrinsic variations of WASP-33 itself. 

Eq.(l) is a generalization of the usual methodology of ra- 
tioing the target star to the sum of the comparison stars. 
Using an equally-weighted sum of the comparison stars 
can be an imperfect divisor, for several reasons. First, the 
comparison stars usually differ in spectral type from the 
target star, and the integral of their flux over our broad 
filter produces a slightly different effective wavelength for 
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each star. The extinction in the stellar signals caused by 
the terrestrial atmosphere can vary strongly with wave- 
length, so different effective wavelengths can degrade the 
correlations between WASP-33 and the comparison stars. 
Second, the comparison stars lie at quite different locations 
on the detector, and higher order effects can be different 
at these different locations. For example, our data exhibit 
a slow drift in the position of all stars over the course of a 
night (about one arc-sec total), caused by differential re- 
fraction between the optical guider and the infrared wave- 
length used for imaging. That small drift, combined with 
small inaccuracies in the flat-field calibration, can produce 
different slow variations for each comparison star. Also, 
slight changes in the degree of defocus can result from in- 
accurate compensation of thermal effects and mechanical 
flexure (Sec. 2.1), and can produce subtle changes in the 
defocused images - hence in the aperture photometry - as 
a function of field position. For these reasons, we believe it 
is good practice to optimally-weight the comparison stars 
using linear regression. 

We compared our optimally-weighted photometry to 
photometry that used an unweighted sum of the compar- 
ison stars. The point-to-point scatter in the two photo- 
metric data sets was about the same, but the baseline was 
much flatter using the optimal technique. Specifically, we 
found that the unweighted sum produced overall slopes of 
0.6% and 0.8% for the nights of 16 and 10 October, respec- 
tively. Using the optimal technique reduced these baseline 
slopes to 0.03% on both nights. 

The linear regression used to solve Eq. (1) can poten- 
tially affect the measured eclipse, because the regression 
attempts to remove all fluctuations in the target star, and 
the eclipse is a fluctuation. In the hypothetical case where 
one of the comparison stars happened to exhibit a noise- 
like decrease at the expected time of the WASP-33 eclipse, 
the regression would overweight that star and the result 
would be to weaken or remove the eclipse. We avoid that 
possibility by solving for the Uj using only out-of-eclipse 
portion of the data (adopting the central phase as 0.5, and 
setting the duration of eclipse to equal the duration of the 
transit.) 

Figure 2 shows the photometry for WASP-33 in the K s 
band on both 10 and 16 October, and the weighted sum 
of the comparison stars (right hand side of Eq. 1) is over- 
plotted as a blue line. The overall rise and fall of these 
signals is due to the airmass dependence of telluric absorp- 
tion, and shorter term fluctuations due to telluric effects 
are also visible. Note also that the weighted sum of the 
comparison stars shows more short-term noise than does 
WASP-33b, because most of the comparison stars arc 1- 
to 2-magnitudcs fainter than WASP-33 (see Table 1). For 
that reason, we smooth the comparison star sum over 10 
points (about 7 minutes in time) before using it to cor- 
rect WASP-33. That degree of smoothing decreases the 
point-to-point noise in the ratio, while still following most 
telluric fluctuations. 

4.2. Warm Spitzer Photometric Analysis 

A dominant effect in Spitzer photometry at both 3.6- 
and 4.5 /im is the presence of intra-pixel sensitivity vari- 
ations. The photometric intensity of a star will depend 

11 http: / /ssc. spitzer. caltcch.edu/warmmission/news/21oct210memo.pdf 



on its position on the detector, and therefore will vary 
with time because of a pointing oscillation in the telescope 
(|Carev et al.ll2010t see lTodorov et al.ll2012l for a recent ex- 
ample). The Spitzer project recently implemented soft- 
ware updates that decrease the amplitude of the pointing 
oscillation, and also decrease its period from 1-hour to 
about 40 minutef^l. This reduces the impact of the intra- 
pixel sensitivity variations on the photometry. Indeed, our 
observations of WASP-33 show the intra-pixel effect to a 
much less degree than many previous observations. 

Many previous investigations have established that the 
intrapixel effect is more dependent on the Y-coordinate 
than on the X- coordinate, and it i s stronger at 3.6- 
than at 4.5 um dKnutson et all 120091: iBeerer et all 120111: 



iDeming et aLll2011at iTodorov et alj|2012ft . We seethe in- 
trapixel effect in our 3.6 /im photometry, and we corrected 
it using a quadratic fit to the photometry as a function of 
the Y-coordinate of the image, with a linear term as a func- 
tion of the X-coordinate (the variation with X is weaker 
than for Y). However, we cannot detect any intra-pixel 
sensitivity variations in our 4.5 /^m photometry. Plots of 
the measured intensity of the star versus both the X- and 
Y-coordinate of the image centroid are essentially scatter 
plots (not illustrated), with no signficant trends. Pear- 
son correlation coefficients of intensity versus X and Y- 
coordinate have values near zero. We also calculated the 
Pearson correlation coefficients for temporal subsets of the 
4.5 /im data, chosen using a moving boxcar window of var- 
ious widths. We can find no significant correlations be- 
tween our 4.5 fjjn photometry and spatial coordinates, dur- 
ing any time period. As an additional check, we repeated 
our photometry using an alternative method to determine 
the position of the star on the detector (center-of- light). 
This method also failed to reveal correlation between po- 
sition and intensity. We therefore use our 4.5 /jm pho- 
tometry for eclipse analysis without applying any spatial 
decorrelation. Note that we do see very small variations in 
the photometry at the 40-minute period corresponding to 
the telescope oscillation (see Sec. 4.4.2). We cannot rule 
out the possibility that the spatial-intensity correlation is 
being obscured by the stellar intensity oscillations. 

In the case of the 3.6 /im photometry, we also used the 
spatia l decorrelation method described by Ba llard et al.l 
pOlOl) . That weighting function method commonly uses a 
time threshold to zero-weight points that lie close in time 
to any given point. When the only intrinsic temporal vari- 
ation is an eclipse or transit, the time threshold is straight 
forward to implement. However, when continuous stellar 
oscillations are part of the desired signal, the weighting 
function time threshold can be problematic. We applied 
a weighting function to the 3.6 /im photometry, using a 
time threshold of zero, i.e. not excluding any points in the 
calculated weights. Fitting to this alternative version of 
the 3.6 /im photometry produced results that differed in- 
significantly from our quoted results (0.34cr, and 0.03cr dif- 
ferences in central phase and eclipse depth, respectively). 
Our quoted results (see Sec. 4.4.2) are based on the poly- 
nominal decorrelation described above, because we believe 
that method is better suited to the nature of these data. 

After decorrelating (or not) the intra-pixel effect, we 
omit some of the initial data for both Spitzer wavelengths. 
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The 3.6 /jm data exhibit an initial transient drift in the 
image position, amounting to about 0.17-pixels, about 4 
times as large as the peak-to-peak variation caused by the 
40-min telescope oscillation. This positional instability is 
accompanied by a similar large transient effect in the pho- 
tometry. We are familiar with the nature of these large 
transient effects based on seeing them in other Spitzer 
data. The relation between intensity and X,Y position 
is not the same during the transient as during the sta- 
ble portion of the time series, because the image traces a 
different region of the pixel. Therefore, we omit the first 
83 minutes of the 3.6 /um data - this being the time for 
the position of the image to stabilize. We see no obvious 
transient effects in the 4.5 /im image position, but the first 
27 minutes of the photometry are anomalously low. As 
a precaution, we omit those initial data from our 4.5 /im 
analysis. 

4.3. Observed Stellar Oscillations 

Figures 3, 4, and 5 show the ground-based WASP-33 
photometry after all correcti ons. WASP-33 is known to 
be an oscillating 6-S c uti star (ICollier-Cameron et al.ll2010t 
iHerrero et all20ll . 15 errero et al.l f|201 lh find a dominant 
oscillation period of 68 minutes, but oscillating stars can 
exhibit different oscillation modes that rise and fall in am- 
plitude. We have therefore decomposed the photome try 
on each night using wavelets (jTorrence fc Compolll998l ) to 
account for possible intra-night changes in the oscillatory 
spectrum. Because the eclipse can significantly affect the 
wavelet analysis, we remove the eclipse prior to computing 
the wavelet power spectra. For the removal, we adopt the 
best-fit eclipse parameters as derived in Sec. 4.4. (This is 
not a circular procedure because the fits for the eclipse pa- 
rameters converge relatively independently of the starting 
estimates for the oscillation periods.) 

Our results show the stellar oscillations very promi- 
nently, but (interestingly) to varying degrees. The J-band 
results are shown in Figure 3; the bottom panel shows 
the corrected photometry and the overlaid blue line shows 
part of the wavelet decomposition - just the lowest order 
components to guide the eye to the general trend. (Note 
that we do not use the blue lines on Figures 3 - 5 in our 
analysis, but we do use the results of the wavelet power 
spectra, shown in the upper panels.) 

The power spectrum in Figure 3 peaks at f46 min- 
utes (frequency 9.9 cy/d), close to twice th e period of the 
68- mi nute oscillation (21.2 cy/d) seen by IHerrero et al.l 
(|2011[ ). No m ention of thi s lowe r frequency mode 
was reported by IHerrero et al.l (|2011l) . but (5-Scuti stars 
are known to exhibit multip le periods of oscillation 
(Bal ona fc Dziembowskil 1 2 llh . wherein the non-radial 
modes c luster around diffe rent radial modes (e.g., n 
1.2..').... iBreger et al.l 120091) ■ Hence, detec t ion o f a dif- 
ferent mode frequency than IHerrero et al.1 f|201 lh is not 
surprising. Figure 4 shows the results in similar format 
for the -RT s -band observations on 10 October. In this case 
a very prominent oscillation is seen with a power spectral 
peak at 71 minutes (20.3 cy/d), in good agreement with 
essentially the same oscillation period (68.56 minute s) that 
was prominent in the results of lHerrero et al.l ([201 ll ) . Note 
also that it can be clearly discerned in the photometry even 
prior to the correction using the comparison stars - see Fig- 



ure 2. The amplitude of this oscillation - about 0.2% in 
flux, i s about twice what was observed by IHerrero et al.l 
(|2011l) in the optical, and we discuss possible reasons for 
this difference in Sec. 5. 

Figure 5 shows the analogous result for the K s band 
observed on 16 October. In this case we see two peaks 
in the power spectrum at 54 and 126 minutes, that are 
probably due to the stellar oscillations as seen on 10 and 
14 October. For each night, and also for the Spitzer data, 
we estimate the total stellar oscillation amplitude from the 
power present in the oscillatory component of the MCMC 
eclipse fit (see Sec. 4.4), and we tabulate those ampli- 
tudes in Table 2. All of these amplitudes are in reasonable 
agreement with the 0.1% oscillation seen in the optical by 
IHerrero et ail (|2011h . except in K s band where the ampli- 
tude is about twice as great, consistently on both 10 and 
16 October. 

Variability in stellar oscillation a mplitudes is sometimes 
obser ved in 5-Scuti stars (e.g., IBreger fc Pamvatnvkhl 
2006), but is not the most likely explanation for our re- 
sults, because of our J-band data. Although J-band shows 
a prominent oscillation near phase 0.82 on Figure 3, quan- 
titative calculation of the average amplitude over that en- 
tire night gives 0.16% (Table 2), only marginally higher 
than the optical value. A close-to-normal oscillation am- 
plitude in J-band, between the two nights of higher am- 
plitude in K s band, seems too much of a coincidence to 
attribute to variability. While we cannot strictly reject 
mode variability, we hypothesize that oscillation ampli- 
tude is associated with wavelength. We discuss possible 
wavelength- related explanations in Sec. 5.2. 

4.4. Eclipse of WASP-SSb 
4.4.1. K s -band Eclipse 

The eclipse of the planet is visible near phase 0.5 in 
the lower panels of Figures 4 and 5. No eclipse is seen 
in Figure 3, as expected for that range of orbital phase. 
The stellar oscillations are a very significant source of con- 
fusion as regards measuring the properties of the eclipse. 
That confusion is mitigated by sufficiently long duration 
of our observed sequences (7.3 hours on 10 October amd 
5.4 hours on 16 October) compared to the periods of oscil- 
lation, and by the fact that the stellar oscillations will not 
be in phase on the two nights, whereas the eclipse should 
repeat at the same orbital phase. In order to exploit the 
latter advantage, we combine the data for 10 and 16 Oc- 
tober into one binned data set, using a bin width of 0.001 
in phase, and we fit models to these binned data. Figure 6 
shows the combined and binned data, with best-fit eclipse 
plus oscillation curves overplotted. 

To obtain the best-fit eclipse depth, we fit a combination 
of eclipse curve plus two sinsusoidal oscillations, using a 
Markov C hain Monte Carlo (MCMC) method, with Gibbs 
sampling ()Fordll2005h . We fit to the combined and binned 
data of Figure 6, in order to gain the advantage of cancel- 
ing the stellar oscillation as much as possible. We gener- 
ate four chains of 10 6 steps each, and we discard the first 
2 x 10 5 steps in each chain. We initialize the MCMC fit 
using sinusoids with periods of 68 and 146 minutes, guided 
by our wavelet power spectra in Figures 3-5. The MCMC 
fit is insensitive to the exact choice of initital periods. As 
long as the initial values contain one period in the ~ 70- 
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minute range, and one period in the ~ 140-minute range, 
the Markov chains rapidly converge to the best-fit values 
for that particular eclipse (this is also true for the Spitzer 
eclipses, Sec. 4.4.2). 

In order to derive accurate MCMC posterior distribu- 
tions, it is essential to re-scale the error bars to yield a 
best reduced x 2 near unity. We find that errors about four 
times larger than the photon noise (for the binned data) 
are necesary to obtain a best-fit reduced \ 2 °f unity. Error 
rescaling factors almost as large as 3 have been required in 
other investigations, ev en for very precise spectrophotom- 
etry (|Bean et al.ll20lil ). An even larger error re-scaling 
in our case is probably a consequence of the complexity 
of the stellar oscillation spectrum, that is only approxi- 
mated using two oscillation periods. In other words, small 
amplitude residual oscillations may masquerade as noise. 

We experimented with adding additional oscillation 
modes to the fit, and indeed this does reduce the ne- 
cessity for error rescaling. However, we can only justify 
initializing the fit using the two periods that we can ob- 
jectively identify in our wavelet power spectra. Other 
modes, if present, cannot be resolved as clear oscillations 
in our ground-based data, instead they appear only as ex- 
tra noise. We also verified that adding a third oscillation 
mode to our MCMC fit does not significantly alter the 
best-fit eclipse depth. 

We vary nine parameters in the MCMC chains: an 
amplitude, period, and phase for two independent sinu- 
soids, an eclipse depth and central phase, and an addi- 
tive constant. We generate the eclipse curve using a new 
version of the iMandel fc Agoll (|2002[) methodology (see 
iDeming et al.ll2011bf ). We adopt system parameters, that 
determine the shape a nd du ration of the eclipse, from 
ICollier-Cameron et al.l (|2010D . excep t for the orb i tal pe - 
riod where we use the value from ISmith et al.l (|201ll ). 
Observations of recent transits in the Czech Exoplanet 
Transit DatabasJ^l sh o w th at the original period from 
ICollier-Cameron et al.l (120101) is too short ; the slightly 
longer period derived by ISmith "eTaD d20ll is more con- 
sistent with recent transit times in the Czech database. 

The top panel of Figure 6 overplots the best-fit curve, in- 
cluding both the eclipse and oscillatory component. The 
lower panel removes the oscillatory component from the 
data, and shows the comparison with the eclipse curve 
alone. Our best-fit eclipse has a depth of 0.0027 ± 0.0002, 
with central phase of 0.4995 ±0.0010. Those errors are im- 
plied by our MCMC posterior distributions. They include 
possible degeneracies between fitted parameters, but there 
can be external sources of error than are not represented in 
the MCMC chains. The most obvious source of such errors 
is the presence of unmodeled stellar oscillation structure, 
as noted above. The most realistic check on the errors is 
to fit each of the two independent nights separately, and 
compare the independent results. Those fits yield eclipse 
depths of 0.0024 ± 0.0002 and 0.0033 ± 0.0004 for 10 and 
16 October, respectively. The best-fit central phases for 
the separate nights are 0.502 ± 0.001 and 0.495 ± 0.001, so 
the difference between two independent nights exceeds the 
formal errors, especially in the case of central phase. We 
adopt the best-fit values determined for the combined and 
binned data, but we assign the errors based on the differ- 

12 http://var2.astro.cz/ETD/ 



ences in the two independent nights. The difference be- 
tween two independent values drawn from a normal error 
distribution is, on average, twice the error associated with 
the average of those two values. So we estimate the er- 
rors appropriate to our average K s -band eclipse depth and 
central phase to be half the difference between the best-fit 
values for the individual nights. Our best-fit eclipse depth 
and central phase is given in Table 3, together with the 
Spitzer eclipse results. 

4.4.2. The Eclipse in Warm Spitzer Bands 

At both Spitzer wavelengths, we solve for the eclipse 
depth and central phase also using an MCMC method, 
with Gibbs sampling. We allow fo r an exponential base- 
line ramp ((Harrington et al.ll2007l ). because we find that 
a linear ramp is inadequate. Inde ed, detailed analy sis of 
very high S/N Spitzer 8/im data (|Agol et al.ll2010l ) indi- 
cates that even a second exponential term can be war- 
ranted. However, we use a single exponential for two rea- 
sons. First, the single exponential ramp alone is more 
comple x than is normally required for these specifi c Spitzer 
bands (|Knutson et all 120091 : iTodorov et alll2012h - expo- 
nential ramps are normally only required at the longest 
wavelength Spitzer bands. Second, our data do not attain 
sufficient S/N to justify the inclusion of a second exponen- 
tial term. 

As in the case of our ground-based data, we calculate 
wavelelet power spectra at both Spitzer bands, and these 
are shown as the upper panels of Figures 7 and 8. Like the 
ground-based data (perhaps a coincidence), we see two sig- 
nificant oscillatory spectral peaks after removing the best- 
fit eclipse. Figure 7 (3.6 fim) shows one of these peaks at 
68 minutes. However, Figure 8 (4.5 /im) shows two peaks, 
39 and 68 minutes. The former is likely to be associated 
with the Spitzer telescope pointing (Sec. 4.2), while the 
latter is clearly due to the stellar oscillation. 

Our first MCMC fits at both Spitzer wavelengths vary 
12 parameters: two sinusoidal terms (each having ampli- 
tude, period and phase), an eclipse of variable depth and 
central phase (i.e., two parameters), the exponential ramp 
(three parameters) , and a zero-point constant in intensity. 
We discard the first 20% of each 10 6 -step Markov chain 
and locate the best fit values from the minimum in \ 2 
(we keep track of \ 2 during the evolution of the chain). 
Subtracting the best fit from the data, we find that the 
residuals have a quasi-sinusoidal shape. These residuals 
represent the unmodeled portion of the data. They can 
be due to red noise in the data created by uncorrected 
instrumental effects, or by imperfect modeling of the stel- 
lar oscillations. Unlike the case of our ground-based data, 
the Spitzer data have sufficient stability and S/N to define 
these small amplitude imperfections in the fit. Both the 
oscillating intensity of the star, and (to a lesser extent) the 
subtlety of red noise caused by Spitzer's improved perfor- 
mance, require some new methodology in their analysis. 
We now introduce that new methodology. 

Following the inital fit using a single MCMC chain of 
10 6 steps, we re-fit by including an additional term to ex- 
plicitly account for the imperfections in the initial fit. We 
calculate residuals for the first MCMC fit by subtracting 
that best fit from the data. We then approximate those 
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residuals using a wavelet decomposition, using N coeffi- 
cients for Morelet wavelets, and we vary N in our sub- 
sequent analysis. For each choice of N, we multiply the 
wavelet decomposition of the residuals by an adjustable 
amplitude, and include that amplitude as a fit parame- 
ter in subsequent MCMC chains. In the limit where N 
equals the number of data points, this procedure would be 
trivial, because the wavelet decomposition would repro- 
duce the residuals exactly, and including those residuals 
in the fit would simply be a fudge, wherein we arbitrarily 
subtract that part of the data not accounted for by the 
model. In the real analysis, this procedure is not trivial. 
With not-too-large N, the wavelet decomposition approx- 
imates only the major features of the residuals, not the 
point-to-point noise. Because the MCMC chain varies the 
amplitude of those major features, correlating the chained 
values of eclipse depth and phase with the amplitude of the 
unmodeled features indicates the degree to which the best- 
fit eclipse depth and phase depend on unmodeled aspects 
of the data. We determine the optimal N by examining 
the noise in the final fits, requiring that it be close to, but 
not less than the photon noise, and have no detectable red 
noise component. This dictates N — 9 for both our 3.6 
and 4.5 /Ltm data. Our final values are based on 4 indepen- 
dent MCMC chains at N = 9, each having 10 6 steps and 
ignoring the first 20% of each chain. 

Using this procedure, we find that the Spitzer eclipse 
depths and central phases are reasonably robust in the 
sense that they do not vary greatly with TV. This is espe- 
cially true for the central phase at 3.6 /im, and the eclipse 
depth at 4.5 /mi. In those cases, the variation in eclipse 
depth and central phase, as we vary N, are consistent with 
the errors implied by the MCMC posterior distributions. 
However, for the converse set of parameters - the 3.6 /im 
eclipse depth and 4.5 /jm central phase - the agreement 
is not as good. In those cases, the differences in best-fit 
values as we vary N are larger than the errors from the 
MCMC posterior distributions, and we adopt the larger er- 
rors consistent with the fluctuations in the best-fit values 
as N varies. 

The derived best-fit eclipse depths, central phases, and 
errors are included in Table 3. We also extract the am- 
plitude of the stellar oscillation in the Spitzer bands (see 
Table 2), based on the total amplitude of the oscillatory 
portion of the MCMC fits, but with the 40-minute por- 
tion subtracted because we attribute that portion to the 
Spitzer observatory. 

4.5. Models for the Atmosphere of WASP-33b 

We interpret our observed planet-star flux contrast of 
WASP-33b using plane-parallel models of the dayside at- 
mosphere of the planet. We use the atmospheric model- 
ing and retr ieval methodology of IMadhusud han fc Seagerl 
(2009], I2010D . The model computes line- by- line radia- 
tive transfer for a plane-parallel atmosphere with the as- 
sumptions of hydrostatic equilibrium and global energy 
balance. The composition and pressure-temperature (P- 
T) profile of the atmosphere are free parameters in the 
model. Since there are as yet insufficient data to con- 
strain abundances in WASP-33b, we adopt solar com- 
position for all elements except car bon, which we allow 
to increase in abundance, following IMadhusudhan et al.l 



(|2011bf) . We explore a variety of temperature profiles, 
especially temperature profiles that are consistent with 
nearly complete absorption of stellar irradiance. The 
models include all the major opacity sources expected 
in hot hydrogen-dominated atmospheres, namely H2O, 
CO, CH4, CO2, TiO, and VO, and collision-induced ab- 
sorption (CIA) d ue to H 2 -H 2 . Our molecular line lists 
are obtained from iFreedman et al.ll2008l iFreedmanl 12 009. 
Rot hman et al.l 120051 iKarkoschka fc Tomaskol 12010 and 
Karkoschkal 201 l| O ur CI A opacitie s are o btained from 
Borvsow et all (11997ft and iBorysowl (|2002f ). A Kurucz 
model ( Castelli k, Kuruczl |2004| ) is used for the stellar 
spectrum, an d the stellar and planetary p arameters are 
adopted from lCollier-Cameron et all (|20 10T) . 

5. RESULTS AND DISCUSSION 

5.1. Stellar Oscillations 

The first and most obvious result from our observations 
is the existence and prominence of the stellar oscillations. 
WASP-33 was already kn own to exhibit oscill ations, but 
the amplitude observed bv lHerrero et al.l (|2011h in the op- 
tical (Johnson R-band) was about 0.001. Our results (Fig- 
ures 3-5, & Figures 7-8) show oscillation amplitudes in 
agreement with the optical, except for K s band where the 
ampitude is about twice the optical value (2.15/Ltm, Ta- 
ble 2), as noted in Sec. 4.3. Because the largest difference 
with the optical amplitude is seen in our ground-based 
(ifs-band) data, we contemplated whether the difference 
could be attributed to errors in our ground-based result. 
An argument against that possibility is the prominence of 
the stellar oscillation in the raw photometry (e.g., upper 
panel of Figure 2). We therefore explore whether proper- 
ties of the stellar atmosphere that may be unique to K s 
band could cause the oscillations to have greater amplitude 
at that wavelength. 

5.2. Stellar Atmospheric Effects 

We here consider the possibility that the larger K s band 
oscillation amplitude as compared to the optical is due to 
the different height of formation for continuum radiation in 
the stellar atmosphere, in concert with height-dependent 
variations in the mode amplitudes. Due to the increase 
in atomic hydrogen bound-free continuous opacity in the 
infrared, our K s band observations of WASP-33 sample a 
greater height in the stellar atmosphere than observed by 
iHerrero et all (|201lh . 

The upward propagation of a pressure-mode oscillation 
in a stellar atmosphere can in principle cause the mode 
amplitude to increase. As a propagating mode encoun- 
ters lower mass density, the wave velocity - and hence the 
temperature perturbation in the compression - increase. 
However, propagation is strongly affected by the strati- 
ficati on of the stellar atmosphere ([Marmolino fc Severinol 
|1991|) . Frequencies less than the acoustic cutoff frequency 
will not propagate, and their velocity amplitude decreases 
with height. To the extent that the temperature ampli- 
tude scales with velocity, it too will decrease with height 
for non-propagating modes. The acoustic cutoff frequency 
is c/2H, where c is sound speed and H is the pressure scale 
height. We calculated the acoustic cutoff frequency and 
other parameters for WASP-33, using a Teff/log(g)/[M/H1 
= 7500/4.5/0.0 model atmosphere from lCastelli k, Kuruczl 
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(2004). We find that the acoustic cutoff corresponds to an 
oscillatory period of about f minute. The much longer 
period oscillations we observe for WASP-33 are therefore 
not propagating, and their amplitudes should decay with 
height. 

The dominant opac ity due to atomic hydrogen 
(jMenzel k Pekerisl [l935h is higher at 2.15 /jm versus the 
optical by a factor of about 1.6, translating to a height dif- 
ference of about 60 km. That is not sufficient to account 
for significant changes in the mode properties, even for 
propagating modes and, as noted above, these modes do 
not propagate. Moreover, if height dependence of the stel- 
lar continuous spectrum were significant to our observed 
amplitudes, then we would expect even larger amplitudes 
in the Spitzer bands than at 2.15 /um, which are not ob- 
served. 

However, there is one unique feature of the K s band 
that may well be responsible for a higher oscillation mode 
amplitude. The Brackett-7 line at 2.165/im is centered 
in the K s bandpass. The strong opacity in that line for 
an A5V star could have a large potential effect on oscilla- 
tion amplitudes. The impact of strong oscillations in the 
line, when diluted over the broad K s band, could poten- 
tially be calculated using techniques beyond the scope of 
this paper (i.e., radiation hydrodynamics). A more direct 
method would be to obtain infared spectroscopy of the 
star and directly measure the oscillations in the infrared 
hydrogen lines. 

5.3. Orbit of WASP-SSb 

Our measured times of central eclipse can in prin- 
ciple determine ec osui for the planet's orbit (e.g., 
iKnutson et al.ll2009l ). Considering the 25 seconds of light 
travel time across the planet's orbit, we expect to find the 
eclipse at phase 0.50024 if the orbit is circular. Weight- 
ing both our ground-based and Spitzer eclipses (Table 2) 
by the inverse of their formal variances, we find an aver- 
age eclipse phase of 0.50044 ± 0.00008, totally dominated 
by the 3.6 /jm eclipse. Thus we find that ecosuj - ap- 
proximated as ir/2 times the phase offset from 0.5 - is 
0.0003±0. 00013. The central phase measurement for these 
eclipses is complicated by the stellar oscillations, so more- 
than-usual caution is needed in the interpretation of the 
measured central phase. Moreover, our value for e cos to 
differs from zero by less than 3er, so our results provide 
little evidence for a non-circular orbit. 

5.4. Atmosphere of WASP-SSb 

Figure 9 shows our result for the eclipse of WASP-33b in 
comparison to two models of its atmosphere, both of which 
agree with the available measurements to date. One of 
these models has a temperature inversion with solar com- 
position, and one has a non-inverte d atmospheric struc- 
ture w ith a carbon-rich composition (|Madhusudhan et al.l 
2011b). Their temperature profiles are shown in the bot- 
tom panel of Figure 9, with the approximate formation 
depths of the four bandpasses overplotted as points. The 
K s bandpass is relatively devoid of strong molecular ab- 
sorption features, and probes the relatively deep planetary 
atmosphere (pressure, P ~ 0.6 bars), relatively indepen- 
dent of the composition of the atmosphere (asterisks on 
Figure 9). (The 3.6 /im bandpass also peaks relatively deep 



in the atmosphere, but has significant contribution from 
higher altitudes, having more molecular absorption than 
does the K s band.) The large eclipse depth we observe in 
K s band (brightness temperature « 3400K) thus indicates 
a hot atmosphere at depth, and a high effective tempera- 
ture for the planet. Both models illustrated on Figure 9 
have hot lower atmospheres (both above 2500 Kelvins) and 
both have inefficient longitudinal energy re-distribution. 
These models are consistent with the observed tendency 
for the most strongly irradiated plan ets to exhibit the leas t 
longitudinal re-distribution of heat (|Cowan k Ago]||20ll . 
We find that these very hot models are nec essary to re- 
produ ce our results as well as the result of ISmith et al.l 
(|2011|). and we concl ude that WASP-33 strengthens the 
ICowan k Ago! (|2011[) result. 

The two models we show are representative of a larger 
set of solutions that explain the data with and without 
therm al inversions. Given that there are 10 model param- 
eters (jMadhusudhan k Seagerll2009L |2010D and only four 
data points, it is not possible to derive a unique model fit 
to the data. We ran large MCMC chains (of ~ 10 6 mod- 
els) with and without thermal inversions, and identified 
regions of co mposition space in each case th at are favored 
by the data (jMadhusudhan k Seagerll2010[ ). 

Both models on Figure 9 are unusual as compared 
for example to the well-o bserved archetype HD 189733b 
(|Charbonneau et a l. 2008) . One model on Figure 9 adopts 
solar composition but with an inverted temperature struc- 
ture (temperature rising with height), while the other 
model has temperature declining with height, but requires 
a carbon-rich composition. We integrate the fluxes of 
each planetary model, and the Kurucz model stellar at- 
mosphere, over the observational bandpasses, and ratio 
those integrals. These band-integrated points are shown as 
squares on Figure 9. The \ 2 values for the models as com- 
pared to_ah|Jojarobserved points (our three measurements, 
plus ISmith et all 1201 ll) are 7.9 for the inverted solar- 
composition model and 2.8 for the non-inverted carbon- 
rich model. The if s -band point at 2.15 /^m favors the non- 
inverted model. Although the difference is not sufficiently 
significant to rule out the inverted model at this time, 
additional eclipse observations in the K s band would be 
helpful to rule out an inverted atmospheric structure. 

In our population of models with thermal inversions, 
several models with slightly different inverted temperature 
structure fit the data almost equally well. However, none 
of them fit the K-band point to within the la errors while 
also fitting the remaining points. The Figure 9 model is 
the best among this set of inverted models. 

In our models without thermal inversions, the best-fit 
model requires a carbon- rich composition (i.e., C/O > 1). 
However, at the 2a level of significance per wavelength 
point, several solar composition models (not illustrated) 
provide an acceptable fit to the data. So although the 
carbon-rich composition is favored, a solar abundance 
composition cannot be absolutely ruled out. 

Our results illustrate the limitations of eclipse photom- 
etry in broad bands, especially for challenging cases like 
planets orbiting oscillating stars. Once we admit the pos- 
sibility of non-solar compositions (because we are largely 
ignorant of true exoplanetary compositions) , the range of 
models that can fit broad-band photometry can be large, 
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in this case extending to inverted and non-inverted mod- 
els with drastically different temperature structure. The 
degeneracy is exacerbated by the relatively small range of 
atmospheric pressures probed by the four bandpasses we 
analyze (points on the bottom panel of Figure 9). Fortu- 
nately, future observations can break th is degeneracy usin g 
HST/WFC3 spectroscopy near 1.4 ^m (jBerta et alj|2012ft . 
The water band near 1.4 /im is sufficiently strong that 
eclipse observations with the Hubble WFC3 grism should 
be feasible. Although water absorption in the carbon-rich 
model is supressed by C/O > 1, the water emission in 
the solar-composition inverted model is predicted to be 
significant, readily detectable near 1.4 /im (see Figure 9). 
Moreover, spectroscopic techniques can potentially probe 
a larger depth range in exoplanetary atmospheres than 
does photometry, because the cores of resolved spectral 
features have strong opacities. Our results therefore illus- 
trate the complementary value of acquiring both broad- 
band and spectroscopic observations of transiting exoplan- 
ets at secondary eclipse. 

6. SUMMARY 

We have analyzed ground-based and Spitzer infrared ob- 
servations of the strongly irradiated exoplanet WASP-33b. 
Our observations span the time of secondary eclipse in K s - 
band and the Warm Spitzer bands at 3.6- and 4.5 /*m. We 
also observed the system for one ground-based night out 
of eclipse in J-band (1.25 /im). Oscillations of the 5-Scuti 
host star are prominent in our data, with a semi-amplitude 
(0.1%) about the same as in the optical. One exception is 
the A^-band, where we find an oscillation amplitude about 
twice as large as seen in the optical. Neither temporal vari- 
ability, nor the variation of stellar continuous opacity with 
wavelength, are likely to account for our A s -band result. 
We speculate that the greater amplitude in A^-band may 
be related to the presence of the Brackett-7 line in the 
bandpass. 

We measure two K s band eclipses, and we base our er- 
rors for the eclipse depth and central phase on the dif- 
ference between the two independent measurements. In 



the case of the Spitzer bands, we measure one eclipse in 
each band. Our Spitzer observations exhibit a relatively 
low level of intra-pixel sensitivity variation as compared to 
previous observations of other exoplanets. However, all of 
our observed eclipses are overlaid by the ubiquitous stel- 
lar oscillations. We adopt an eclipse model comprised of 
the eclipse shape, plus two sinusoids to account for the 
stellar oscillations. We fit the model to the photometry 
using an MCMC method. The Spitzer photometry is of 
sufficient quality that we are able to implement a wavelet- 
based technique to define fluctuations in the data that are 
not accounted for by our fitting procedure. Including that 
structure as a parameter in subsequent MCMC fits, we 
thereby explore how the eclipse depth and central phase 
vary as a function of the amplitude of the unmodeled struc- 
ture. The derived eclipse depth and central phase of the 
Spitzer eclipses are not strongly sensitive to unmodeled 
structure in the light curves, but this procedure does al- 
low a realistic evaluation of the errors. 

Using the million- model ap proach pioneered by 
iMadhusudhan Ez Seager! (|2009ft andlMadhusu dhan k, Seageil 

we explore the range of atmospheric temperature 
structures and compositions that are consistent with our 
eclipse observation s, plus the 0.9 /im eclipse observed by 
Smi th et al.l (|201l[) . We find two possible atmospheric 
models. One has an inverted temperature structure, and 
solar composition, and the second has a non-inverted 
temperature structure with a carbon-rich composition. 
We point out the value of infrared eclipse spectroscopy 
using moderate resolving power to detect (for example) 
the water band near 1.4 /im. Spectroscopic detection of 
water emission or absorption at eclipse would break the 
degeneracy between the two possible models. Although 
the temperature structure of WASP-33b is currently am- 
biguous, both models re-emit a large fraction of incident 
stellar irradi ation from their day s ides, strengthening the 
hypothesis of lCowan k Agol! (| 2 1 lh that the most strongly 
irradiated planets circulate energy to their night sides with 
low efficiency. 
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Fig. 1. — Results of our dome-test observations and procedure to correct for the instrument-related signal instabilities. Top panel: original 
data (red points, connected with line) and corrected data (black points) for our test observations at the WASP-33 position. Middle panel: 
Derivative of the original data from the top panel, with the red points marking outliers in the derivative that are corrected by our methodology. 
Bottom panel: Ratio of the corrected signal (black points) at the WASP-33 position to the sum (not illustrated) of the corrected signals at 
the comparison star positions. Note the greatly expanded intensity scale on the lower panel. These test observations were made observing 
K-band background with the telescope mirror and dome both closed. The cadence of these 2-second exposures was about one per minute, 
with periodic telescope motions in 0.5-hour increments of right ascension from -4 hours to +2 hours, at constant +32-degree declination. The 
instabilities (e.g., at exposure 38) always correspond to times of telescope motion, but not every telescope motion causes an instability. 
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Fig. 2. — Photometry of WASP-33 (black points) in the its-band on 10 and 16 October, 2011, spanning the time of secondary eclipse. 
These data have been corrected for instrumental instabilities, but not normalized using the comparison stars. The blue line is the weighted 
sum of the comparison stars, i.e., the right hand side of Eq. (1), smoothed over 7 points (about 5 minutes of time). 
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Fig. 3. — Bottom panel: Observations of WASP-33 in the J-band on 14 October, 2011, when no transit or eclipse occurred (near phase 
0.8). These data have been corrected for instrument instabilities and also corrected using the comparison stars. Note the oscillatory behavior 
of the star, with an amplitude of about 0.16%. The blue line is a de-noised version using wavelets; it only incorporates the first few wavelet 
coefficients so as to guide the eye to the general variations. Top panel: Wavelet power spectrum of the data from the bottom panel. The 
oscillatory power peaks near 146 minutes. 
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Fig. 4. — Bottom panel: Observations of WASP-33 in the i^s-band on 10 October, 2011, spanning a secondary eclipse. These data have 
been corrected for instrumental instabilities and also corrected using the comparison stars. The blue line is a de-noised version using wavelets; 
it only incorporates the first few wavelet coefficients so as to guide the eye to the general variations. Note the very prominent stellar oscillation 
with period near 71 minutes. Top panel: Wavelet power spectrum of the data from the bottom panel, showing the oscillation peak power at 
71 minutes. 
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Fig. 5. — Bottom panel: Observations of WASP-33 in the i^s-band on 16 October, 2011, spanning a secondary eclipse. These data have 
been corrected for instrumental instabilities and also corrected using the comparison stars. The blue line is a de-noised version using wavelets; 
it only incorporates the first few wavelet coefficients so as to guide the eye to the general variations. Top panel: Wavelet power spectrum of 
the data from the bottom panel, showing oscillatory power peaks at 54 and 126 minutes. 
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Fig. 6. — Top panel: Eclipse of WASP-33b based on combining the if s -band observations from 10 and 16 October, and binning to a phase 
resolution of 0.001. The blue line shows the result of fitting to the eclipse, and the sum of two oscillation modes, via Markov chains. Bottom 
panel: Data from the top panel with the oscillatory portion of the fit removed and compared with the best- fit eclipse (red line). The scatter 
per binned point is about 0.0012, indicated by the inset point with error bars. The best-fit eclipse depth is 0.0027 ± 0.0002, but comparison 
of the two individual nights (not illustrated) indicates a greater error in the depth (±0.0004, see Table 3). 
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Fig. 7. — Middle panel: Spitzer observations of WASP-33 at 3.6/im, spanning a secondary eclipse. The overplotted blue line is the best-fit 
solution from our MCMC analysis, including structure defined by our wavelet analysis (see text). The dashed blue line omits the wavelet- 
defined structure and uses only the pure oscillatory portion, plus eclipse. Bottom panel: Data from the middle panel with the oscillatory plus 
wavelet portion removed, showing the best-fit eclipse curve (red line). Top panel: Wavelet power spectrum of the data points from the middle 
panel. The peak near 68 minutes is due to oscillations of the star. 
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Fig. 8. — Middle panel: Spitzer observations of WASP-33 at 4.5 fim, spanning a secondary eclipse. The overplotted blue line is the best-fit 
solution from our MCMC analysis, including structure defined by our wavelet analysis (see text). The dashed blue line omits the wavelet- 
defined structure and uses only the pure oscillatory portion, plus eclipse. Bottom panel: Data from the middle panel with the oscillatory plus 
wavelet portion removed, showing the best- fit eclipse curve (red line). Top panel: Wavelet power spectrum of the data from the middle panel; 
the main peak near 68 minutes is due to oscillations of the star, and a secondary peak near 39 minutes is due to a pointing oscillation in the 
telescope. 
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Fig . 9. — Top panel: Comparison of our measured eclipse depths (red-orange points at 2.15, 3.6, and 4.5 fim, and including ISmith et al 
201ll at 0.9 fim) with two models for the atmosphere of WASP-33b: the green line is a carbon-rich non-inverted model, and the violet line is 
a solar composition model with an inverted temperature structure. The squares are the values expected by integrating the planetary fluxes 
over the bandpasses. Bottom panel: Pressure-temperature profiles for the models whose emergent spectra are shown in the top panel. The 
peaks of the contribution functions for each bandpass are plotted as points on the pressure-temperature profiles. Squares are 0.9 fim, asterisks 
are 2.1 fim, triangles are 3.6 fim, and diamonds are 4.5 fim. 



Eclipse of WASP-33b 



21 



Table 1 

Comparison stars used for J- and JG-band photometry 



i 


2MASS Designation 


J-mag 


K-band mag 


1 


02263012- 


-3724227 


7.416 


6.811 


2 


02273482- 


-3728084 


7.886 


7.380 


3 


02263566- 


-3735479 


8.541 


8.304 


4 


02262831- 


-3732264 


9.376 


8.579 


5 


02265167- 


-3732133 


9.958 


9.937 


6 


02271232- 


-3728286 


9.808 


9.536 



Note: All stars were used for J-band photometry, but for if s -band we omitted star 1 on 10 October because it was too 
bright, and we omitted star 5 on 16 October because it was not sufficiently above the higher thermal background on that 
night. Also, star 6 was too faint for if s -band on both 10 and 16 October. For reference, WASP-33 has J-mag=7.581 and 
K-mag=7.468. 
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Table 2 

wasp-33 stellar infrared oscillation amplitudes and time scales 



UT Date (2011) 


Wavelength (/xm) 


Oscillation amplitude 


Time scale(s) 


Observatory 


26 March 


3.6 


0.0013 


68 min 


Warm Spitzer; Figure 8 


30 March 


4.5 


0.0013 


68 min 


Warm Spitzer; Figure 7 


10 October 


2.15 


0.0023 


71 min 


Kitt Peak 2-meter; Figure 4 


14 October 


1.25 


0.0016 


146 min 


Kitt Peak 2-meter; Figure 3 


16 October 


2.15 


0.0021 


54, 126 min 


Kitt Peak 2-meter; Figure 5 



Note: 



Smith et al 



(|2011l ) found oscillation periods between 42 and 77 minutes, and 



Herrero et al 



(|2011l ) found a dominant period near 68 minutes. 
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Table 3 

wasp-33b secondary eclipse depths, time of central eclipse, and orbital phase 



Wavelength (fim) Eclipse Depth Time as BJD(TDB) 

2.15 0.0027 ±0.0004 2455844.8156 ± 0.0040 

3.6 0.0026 ±0.0005 2455647.1978 ± 0.0001 

4.5 0.0041 ± 0.0002 2455650.8584 ± 0.0005 



Phase Comment 

0.4995 ± 0.0035 average of 10 and 16 October; Figure 6 

0.50041 ± 0.00008 Warm Spitzer; Figure 7 

0.5012 ± 0.0004 Warm Spitzer; Figure 8 



Note: For the ground-based eclipse at 2.15 /im the barycentric time of central eclipse is based on the average phase for both 10 and 16 October, then converted 
to the central BJD value for the 10 October eclipse. 



